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ABSTRACT 

We revisit the analysis of the Non- linear Thin Shell Instability (NTSI) numerically, including magnetic 
fields. The magnetic tension force is expected to work against the main driver of the NTSI - namely 
transverse momentum transport. However, depending on the field strength and orientation, the 
instability may grow. For fields aligned with the inflow, we find that the NTSI is suppressed only 
when the Alfven speed surpasses the (supersonic) velocities generated along the collision interface. 
Even for fields perpendicular to the inflow, which are the most effective at preventing the NTSI 
from developing, internal structures form within the expanding slab interface, probably leading to 
fragmentation in the presence of self-gravity or thermal instabilities. High Reynolds numbers result 
in local turbulence within the perturbed slab, which in turn triggers reconnection and dissipation of 
the excess magnetic flux. We find that when the magnetic field is initially aligned with the flow, there 
exists a (weak) correlation between field strength and gas density. However, for transverse fields, this 
correlation essentially vanishes. In light of these results, our general conclusion is that instabilities are 
unlikely to be erased unless the magnetic energy in clouds is much larger than the turbulent energy. 
Finally, while our study is motivated by the scenario of molecular cloud formation in colliding flows, 
our results span a larger range of applicability, from supernovae shells to colliding stellar winds. 
Subject headings: instabilities — MHD — turbulence — methods: numerical — ISMxlouds — ISM: 
magnetic fields 



1. MOTIVATION 

Shocks and shells exist abundantly in the interstel- 
lar medium (ISM). Driven by supernovae, expanding 
HH-regions, gravitational flows or wholesale cloud col- 
lisions, they not only strongly influence the ISM dy- 
namics, but also affect its chemistry. However, struc- 
tural analyses of the ISM show spectral indices closer 
to t he Kolmogorov value of incompressible turbulence 
(see lElmesrreen fc Scale! 120041 for a review) . Indepen- 
dently of the problem of how well these indices con- 
strain the type of turbulence, there are plenty of physical 
mechanisms to explain the closeness to the Kolmogorov 
value, rangin g from the intrinsic nature of MHD tur- 
bulence (e.g. iGoldreich fc Sridharl|l995t iBoldvrev et al.1 
2002, iCho fc Lazarianl 12003) to the conversion of com- 
pressible to solenoidal modes (jFalgarone et alJ 11994 
lElmeereen fc Scaldl2l)04|) . 

It is the latter mechanism that motivated this study. 
In the absence of shear flows (oblique shocks) and ther- 
mal instabilitie s, the Non-linear Thin Shell Instability 
l)Vishniacll 994n provides a natural mechanism to convert 
compressible motions into solenoidal ones. The NTSI is 
a rippling instability, relying on transverse momentum 
transport due to bends in the collision interface of two 
opposing flows. It is likely to arise in a wide range of envi- 
ronments, from colliding stellar winds, supernova shells, 
colliding HI streams/clouds, to galaxy mergers. 
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The NTSI has b een widely studied numerically (see 
iHeitsch et akl 120061 for a summary of the literature), 
mostly focusing on the effects of self-gravity and ther- 
mal instabilities. Our interest in the NTSI comes 
from the role it plays in the e volution of molecu- 
lar cl ouds in colliding HI flows <|Heitsch et all 120051 
2006; Ivlzo' uez-Semadeni et akl 120061) . but the results 
have a wider applic ability. With the exception of 
iKlein fc Woods! l)1998j) . who included the magnetic pres- 
sure term in their study of cloud collisions, work on 
the NTSI has so far neglected the effect of magnetic 
fields. However, fields could have a deciding influence 
on the evolution of a sh ock-bounded slab. Motivated by 
the numerical models oflVazc iuez-Semadeni et al.l (1995) 
and iPa ssot et akj lll995| ). lHartmann et alJ l)200lj) and 
iBergin et al.1 l)2004j) suggested that fields could, in fact, 
lead to a selection effect for molecular cloud formation: 
clouds can only form if the fields are aligned with the 
flows assembling the gas. 

As a first ste p, we revisit the isothermal analysis of 
Vishniac (1994) numerically, and study the evolution of 
the NTSI in a two-dimensional, magnetohydrodynami- 
cal environment. While the general expectation ( W2 2.1f> 
is met in the laminar case, namely that magnetic fields 
can efficiently damp the NTSI, the geometry of the rip- 
pled interface induces non- ideal MHD effects, requiring 
a numerical method capable of handling such effects in a 
stable and accurate way ( ^2 2.2J) . However, we find that 
the exact amount of dampening crucially depends on the 
field orientation and strength (Q. This is especially true 
in the turbulent case, where turbulent reconnection in- 
side the over-pressured slab leads to a pressure deficit, 
thus compressing the gas even further. Finally, we show 
that the correlation between field strength and gas den- 
sity is at best weak, even in the case of fields perpendic- 
ularly oriented with respect to the inflow, for which the 
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field would be expected to scale linearly with the density. 

2. PHYSICS AND NUMERICS 

2.1. Physics 

The growth rate of the NTSI is mostly controlled by krj, 
the product of the wave number of the slab perturbation 
k, and the amplitude of the slab's initial displacement 
r\ (equivalently, the amplitude of the collision interface's 
geometrical perturbation). The instability is driven by 
lateral transport of longitudinal momentum, i.e. if the 
inflow is parallel to the x direction, and the slab is in the 
j/-z-plane, x-momentum is transported laterally in y (and 
z), collecting at the focal points of the perturbed slab. 
The efficiency of lateral momentum transport is key to 
the development of the instability, since it is the imbal- 
ance of ram pressure at the focal points that eventually 
propels matte r forward , drivin g the growth of the slab's 
perturbation. iVishniacI (|1994fl derived a growth rate of 

u^cMkrj) 1 / 2 , (1) 

with the sound speed c s . iBlondin fe Marks! l)1996j) found 
that at constant r\ and for small fcs, equation yields 
only a lower limit, while for large fcs, the analytical 
growth rates agree well with the numerical results. The 
reason for this seems to lie in the efficiency of deflecting 
the incoming flow: for small fcs, a small fraction of the 
incoming flow's momentum is converted to lateral mo- 
tions, while a large part compresses the slab (depending 
on the equation of state, this could lead to an increase 
in energy losses). 

After the initial growth-phase (eq. PP), the NTSI 
reaches saturation through two mechanisms: (i) expan- 
sion of the slab which stops the lateral momentum trans- 
port by preventing the inflow from reaching the focal 
points, and (ii) shear flow (Kelvin-Helmholtz) instabil- 
ities (KHI) in regions of the slab connecting the focal 
points. The KHIs both generate inner structure and re- 
vive the slab's expansion. Note that strong cooling (not 
modeled he re) can also suppre ss the NTSI via early frag- 
mentation l|Hueckstaedtll2003[K 

Qualitatively, we expect magnetic fields to prevent the 
NTSI and subsequent KHI-modes from occurring. How- 
ever, the detailed quantitative extent of the damping 
should depend on the orientation of the field with respect 
to the inflow. Indeed, fields aligned with the inflow resist 
instabilities via the magnetic tension force, and therefore 
should be more efficient in suppressing the NTSI when 
krj is small, even though the strong pairwise field rever- 
sals arising from the opposed shear velocities along the 
slab (see W3 3.1(1 could trigger reconnection. On the other 
hand, fields perpendicular to the inflow (but still in the 
2D plane), primarily prevent instabilities from growing 
because of the magnetic pressure term in the Lorentz 
force, and to a lesser extent, because of magnetic tension. 
For the sake of completeness, we mention that the third 
possible field configuration in 2D, i.e. the one in which 
the magnetic field is perpendicular to the flow plane, is 
irrelevant for this study since in that case the gas behaves 
as a system with an adiabat ic exponent of 2 for which 
the NTSI cannot be excited l)Vishniadll99^ . 

2.2. Numerics 



The magnetohydrodynamical scheme 5 is based on 
a conservativ e gas- kinet ic flux-splitt i ng m ethod, intro- 
duced by Xu (1999) and lTang fc Xul I200T1) and derived 
from the lst-order BGK ijBhatnagar et al.lll954j) model. 
Representing the velocity distributions as Maxwellians 
in each cell, fluxes across cell walls are derived from the 
differences in the velocity moments of Maxwellian distri- 
butions reconstructed at the cell walls. The reconstruc- 
tion is second order in space using MUSCL limiters, and 
it allows a fast and consistent way to implement vis- 
cosity and Ohmic resistivi ty in the form of dissipative 
fluxes {Heitsch et alJ 12004^ at close to zero extra com- 
puting cost, while preserving the time order of PRO- 
TEUS since the dissipative terms are not simply added as 
source terms but are part of the flux computation. This 
allows us to control dissipation in a physical manner, 
without having to rely on numerical dissipation to ter- 
minate the turbulent cascade at grid scale. Total energy 
is conserved at machine-accuracy level for an adiabatic 
equation of state. In the isothermal version which we 
are using in the total energy equation is not evolved. 
PROTE US uses a 2nd orde r TVD Runge-Kutta time 
stepping l|Shu fc Osherf ll988) for the MHD equations to 
achieve 2nd order temporal accuracy. Fluxes are up- 
dated in time-unsplit fashion, i.e. flux updates for spa- 
tial directions are computed using the initial conditions 
of the current time step. In order to ke ep V • B = 0, 
PROTEUS e mploys a Hodge projection {Balsaralll99l 
iZacharv et allfl99# . The code is fully message passing 
interface (MPI) parallelized. 

With PROTEUS, one may switch between the MHD- 
solver previously described and a purely hydrodynamical 
solver based on the 2nd-order BGK model. The latter 
implementatio n has been introduce d and extensively dis- 
cussed bvlPrendergast fc Xul (ll993T).lS lvz fc P rendergastl 
ljl999j) . lHe"itsch et all l|2006D and lSlvz et all (|2006l) . and 
we therefore refer the reader interested in implementa- 
tion details to these papers. 

One-dimensional shock tests and t he Orszag-Tang vor - 
tex have already been discussed bv iTang fc XiJ {2000), 
hence in what follows, we focus on three other MHD test 
cases. The two first ones, i.e. the propagation of a linear 
Alfven wave under resistive damping f H2 2.2 2.2. ljl and 
the current sheet evolution ( H2 2.2 2.2.2J) are both meant 
to test the resistive fluxes, while the third one, i.e. the 
advection of a field loop ( W2 2.2 2.2.3fl is a geometry test. 
A detailed study of the magnetized Kelvin -Helmholtz In- 
stability is under way (Palotti et al. 2006). 

2.2.1. Propagation of a linear Alfven Wave 

This one-dimensional test checks the resistive flux 
implementation as well as the accuracy of the overall 
scheme. A linear Alfven wave under weak Ohmic dissi- 
pation is damped at a rate of 

Ui = \^k 2 , (2) 

5 We called the scheme PROTEUS, under which name we will 
refer to it subsequently. Proteus is a lesser god in Greek mythology, 
also known as "The Old Man of the Sea" . He lives in the sea off the 
coast of Egypt and can see things in the past, presence and future, 
but is very unwilling to share his knowledge. In order to evade 
questions, he has the ability to change his appearance. However, if 
you manage to catch and hold him, he will assume his true shape 
and answer your questions. 
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Fig. 1. — Damping rate (eq. [2]) of a linear Alfven wave against 
Ohmic resistivity for re = 1, 2, 4. The resolution is N = 64. Errors 
of the measured damping rates are smaller than the symbol sizes. 
Lines denote the analytical solution. 



where An is the Ohmic resistivity, and k = 2t:k/L is 
the wave number of the Alfven wave, with k £ N. The 
strongly damped case, where the decay dominates the 
time evolution, is physically uninteresting for our appli- 
cation, since the Ohmic resistivity is mainly used to con- 
trol numerical dissipation. Figure ^ shows the damping 
rate against Ohmic resistivity An for re = 1, 2,4 at a grid 
resolution of N = 64. 

From Figure ^ it is clear that, as one diminishes the 
value of An, there comes a point when the numerical re- 
sistivity of the scheme becomes comparable to the phys- 
ical one, causing the measured damping rate to flatten 
out and depart from the analytical solution. For n — 4 
and An =0.1, the wave decays too quickly to allow a re- 
liable measurement, and the system enters the strongly 
damped branch of the dispersion relation. However, we 
emphasize that, even for this high value of k in light 
of the modest resolution we used, the resistivity range 
available to PROTEUS spans three orders of magnitude. 

2.2.2. Current Sheet 



This test is taken from Haw lev fe Stond (^995) and 
the ATHENA test suite 6 . A square domain of extent 
—0.5 < x, y < 0.5 and of constant density no = 1 and 
pressure po is permeated by a magnetic field along the 
y direction such that i?j,(|a;| < 0.25) = y/An, and B y — 
—^/Ak elsewhere. This results in two magnetic null lines, 
which then are perturbed by velocities v x = Asm(2Tty). 
The goal is to find the pressure po and velocity amplitude 
A for which the code crashes. The main problem is - 
especially in conservative schemes - that the resistive 
decay of the field leads to strong localized heating that 
in turn generates strong magnetosonic waves. Thus, the 
smaller po and/or the larger A, the harder the test. We 
chose po = 0.1 and A — 0.3 fairly close to the "standard" 
values quoted on the ATHENA web site 6 , po — 0.05 and 
A = 0.1. Here, we use an adiabatic exponent of 7 = 5/3 
and employ the conservative formulation of the scheme. 

6 http: //www. astro .princeton. edu/ 
~j stone/tests/f ield-loop/Field-loop .html 
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Fig. 2. — Total magnetic energy against time for the current 
sheet test. A finite resistivity Aq helps stabilize the code. 



The test turns out to be very hard for PROTEUS, 
because of its low numerical diffusivity Figure [5] sum- 
marizes the test results in the form of the total magnetic 
energy against time. The energy evolution splits into 
two branches: one corresponding to the lower-resolution 
models at all resistivities, and the other representing the 
higher-resolution models at An = and An = 10~ 5 . At 
An = 10~ 4 , the field decay has converged: both resolu- 
tions give the same curve. At N = 256 2 and An = 0, the 
code crashes around t = 6. All other models run up to 
t = 10 and further. 

The result confirms the discussion by iTang fc Xul 
(2000), namely that while the conservative gas-kinetic 
flux splitting method in the BGK-formalism performs 
well for high-/? plasmas (with (3 defined as the ratio of 
thermal over magnetic pressure), it might not be the 
method of choice for low-/3 plasmas, i.e. magnetically 
dominated systems. Note, however, that this is mainly 
a consequence of the scheme's conservative formulation. 
We experimented with a non-conservative version (i.e. 
just evolving the internal energy instead of the total en- 
ergy), which was stable for lower po and higher A, albeit 
at the cost of a less accurate total energy evolution. 

2.2.3. Advection of a Field Loop 

A cylindrical current distribution (i.e. a field 
loop) is advected diagonally across the simulation do- 
main. We follow the implementation presented by 
iGardiner fc Stond (|2005f) and t he ATHENA test suite 6 , 
based on an earlier version bv iToth fc Odstrcill {1996). 
Density and pressure are both initially uniform at no = 1 
and p = 1, and the fluid is described as an ideal gas 
with an adiabatic exponent of 7 = 5/3. The square grid 
ranges from —0.5 < x < 0.5, and the loop is advected at 
an angle of 30 degrees with respect to the x-axis. Thus, 
two round trips in x correspond to one crossing in y. The 
amplitude of the field loop is set to values 10~ 3 , 10 -2 and 
10 _1 , with an initial radius of i?o = 0.3. Figure shows 
the initial current distribution with the magnetic field 
vectors over-plotted {left), and the final current distribu- 
tion after two time-units measured in horizontal crossing 
times. The overall shape is preserved, although some ar- 
tifacts are visible at the upper rim of the loop. These 
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Fig. 3. — Current density for field loop advection test (see 
t|2 2,2 2.2.31 . Left: Initial condition. Bight: After two horizontal 
crossings. The grid resolution is 256 2 . 
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Fig. 4. — Normalized magnetic energy against time (in units 
of horizontal crossing time), for two resolutions and three field 
strengths. The highest amplitude leads to waves when perturbed 
by advection. At a resolution of 128 2 , the magnetic energy has 
decayed by 3.3% after two crossing times. The fit decay times r 
are indicated for each model. 



results concerning the shape are qualitatively similar to 
those posted on the above mentioned website. Note that 
we show the current density, since the artifacts do not 
show up in the magnetic energy. This test uses An = 0. 

Figure ^ presents a quantitative diagnostic of the be- 
havior of the code as it tracks the magnetic energy de- 
cay against the simulation time, in units of horizon- 
tal crossing times. The initial energy is normalized to 
one. Line styles denote different amplitudes of the field 
strength. The solid lines correspond to the case given by 
iGardiner fc Stone) l|2005D . the dashed and dash-dotted 
lines denote cases with larger amplitudes. At an ampli- 
tude of A — 10~ 3 and a resolution of 128 2 , the magnetic 
energy decays by 3.3% over two horizontal crossing times. 
The ATHENA website 6 quotes a decay of 3.5% with a 
256x148 grid, using a Roe solver and 3rd order recon- 
struction. 

In summary, these numerical test cases demonstrate 
that PROTEUS models dissipative MHD effects accu- 
rately, due to a low intrinsic numerical diffusivity that 
compares well with that of higher-order Godunov meth- 
ods. Furthermore, it can advect geometrically complex 
magnetic field patterns properly, and is well suited in its 
energy conserving form to model MHD flows with (3 > 1. 



2.3. Initial Conditions 

To remain as close as possible to iVishniacl l)1994j) , we 
will use the isothermal version of PROTEUS in the fol- 
lowing. The initial conditions are similar to those dis- 
cussed in lHeitsch et alJ 1)2005112006]) . Two uniform, iden- 
tical flows in the x-y computational plane initially collide 
head-on at a sinusoidal interface with given wave num- 
ber k y and amplitude rj. The field is either aligned or 
perpendicular to the inflow, but in both cases in the x- 
y plane. For the standard runs, we used a rectangular 
grid with an extent of 88 pc in x and 44 pc in y. Field 
strength as well as viscosity and resistivity are varied. 
The grid resolution varies between N x x N y — 256 x 128 
and 2048 x 1024 by factors of 2 in linear resolution. The 
isothermal sound speed is cs = 5.3 km s _1 , the Mach 
number of the instreaming gas is M. = 4, and the inflow 
density is set to no = 1 cm -3 . Thus, in the code unit 
system, the Alfven speed in the inflow region is given by 
ca = B, the magnetic field strength. 

3. RESULTS 

We give a rough estimate for the field strength required 
to prevent the excitation of the NTSI in t|3 3.11 The mor- 
phology of the instability naturally depends strongly on 
the field orientation. We present some examples in H3 3.21 
Because of the strong shear flows, the explicit control of 
dissipation is crucial in reaching numerical convergence. 
This can be further quantified by monitoring the growth 
rates f fl3 3.3(1 . Finally we show that the geometry of 
the flow and magnetic field strongly influence the field- 
density relation f H3 3.4(1 . 

3.1. Estimate of Threshold Field Strength 

A very rough estimate of the threshold field strength 
preventing the excitation of the NTSI can be derived by 
simple pressure considerations. Figure \5\ gives a sketch 
of the simplified situation. Only one half of the slab in 
the vertical direction is shown. The slab is displaced by 
?; in the horizontal direction around the (dotted) center 
line. The angle between the slab and the symmetry line 
measured at point is given by a, with tana rts 2rjk/-K. 
Gas is streaming in horizontally from the left and the 
right with velocity u, and the magnetic field B is aligned 
with the inflow, pointing to the right. Incoming flow with 
positive velocities exerts a pressure at point 0, which can 
be split into a normal component 0D and a tangential 
component OF = 0E sin a. The tangential component 
corresponds to the ram pressure exerted by material de- 
flected by the slab and sliding along OF. To a zeroth 
order approximation, the component of this ram pres- 
sure perpendicular to B is available for bending the field 
lines, thus 

pu 2 sin a cos a ps B 2 /2. (3) 

The maximum is reached for a = 7r/4, in which case 
pu 2 w B 2 . For an inflow velocity of |u| = 4c s , we thus 
require a field strength corresponding to an Alfven speed 
of ca ~ 4c s to suppress the NTSI. 

3.2. Morphology 

We begin with our standard runs, in order to com- 
pare to Vishniac's analysis. These are the "laminar" 
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Fig. 5. — Sketch of force geometry for estimating the threshold 
field strength. A thick solid line marks the interface of the colliding 
flows (the "slab"). The magnetic field is represented by horizontal 
long short-dashed arrows, and inflow velocities by short solid ones. 
Other (force) vectors are defined in the text. 



cases ( , M3 3.2 3.2~l"|) . i.e. cases that do not develop tur- 
bulent substructure. The turbulent case and the rele- 
vance of fixed physical dissipation scales are discussed in 

3.2.1. Laminar Case 

The left column of Figure displays a model sequence 
in resolution for the hydro runs, corresponding to Vish- 
niac's analysis. The top panel shows the initial condition 
(strictly, just after t = 0), The increasing overall ampli- 
tude of the slab signifies that the NTSI is clearly at work. 
Most of the gas is collected at the focal points, and by 
the end of the simulations, the system is close to satu- 
ration. The lowest resolution run differs from all others 
in that it is the only one for which numerical conver- 
gence has not been achieved. The two highest resolution 
runs (N x = 1024,2048, lower two panels), on the other 
hand have converged even in detail. Note that the vis- 
cosity and resistivity provide fixed physical dissipative 
scales, independent of the resolution. Thus, the model 
at N x — 256 is not resolved with respect to these dissi- 
pative scales, while the high-resolution models are. 

This is also true for the magnetic runs (right column 
of Figure |BJ) where the field has slowed down the growth 
of the NTSI. The resulting slab is also more structured 
than in the pure hydro case. High-density regions, es- 
pecially thin filaments, coincide with regions of field re- 
versals (loss of magnetic pressure support). To see this, 
one can compare the gas density to the magnetic energy 
(Figure 0) maps. The center column corresponds to the 
right column of Figure EJ Not only do the field rever- 




FlG. 6. — Left column: Logarithmic density maps of hydrody- 
namical models at t = 3.8 Myr corresponding to the end of simula- 
tion, resolution increasing from top to bottom (N x = 256 to 2048 
by factors of 2). The two highest resolutions have converged even 
in detail. Right column: Same as left, but for the magnetic models, 
where the field is aligned with inflow and ca/c s = 1.0. Again, the 
two highest resolutions have converged. 



sals lead to dense structures, but magnetic waves arise. 
The left column of Figure [7| shows the same resolution 
sequence at half the field strength, i.e. ca/c s = 0.5, 
while the right column stands for ca/c s — 2.0. Higher 
field strengths not only reduce the growth rate of the 
NTSI, but also suppress the internal turbulent structure 
visible for ca/c s = 0.5. Numerical convergence is more 
easily achieved with higher field strength, which is an- 
other indicator that turbulence plays a minor role, i.e. 
we remain safely entrenched within the laminar regime. 
This is slightly different in the weak-field case where the 
two highest resolution runs have only mildly converged. 
Here, the field starts to get too weak to prevent the ex- 
citation of KHI modes. 

A variation on the theme is shown in Figures [S] and [§] 
where the field is oriented vertically this time around. In 
a ID geometry, the magnetic pressure would prevent the 
gas from efficiently acc umulating to form high-density 
regions (i.e. clouds, see He rein et al.ll2004|) . The system 
would behave as if the gas had an adiabatic exponent of 
2. This is the situation traced out by the models shown in 
the right column of Figures |HI and because of the stiff- 
ened equation of state, high (flux-)density regions expand 
faster. However, high density is found at the focal points, 
hence the initial perturbation of the slab is smoothed out, 
and the slab ends up with plane-parallel shock fronts 
(due to the fields, there are some waves inside the slab, 
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Fig. 7. — Logarithm of magnetic energy at t = 3.8 Myr. Left: Four resolutions increasing from top to bottom (N x = 256 to 2048 by 
factors of 2) at ca/c s = 0.5. Center: For c^/c a = 1.0. Right: For ca/c s = 2.0. 



though). Reducing the field strength (left column) by 
a factor of 2 however changes the situation drastically. 
Although the NTSI is only weak, high-density filaments 
start to form in the thick slab, again at the locations 
of field reversals. Thus, a transverse field is much less 
of an inhibiting factor for substructure or high-density 
region generation than what would be expected from a 
ID-argument. For ca/c s = 1.0, the fast magnetosonic 
modes have already reached the boundaries, causing the 
geometric patterns visible in the density plot. 

3.2.2. Turbulent Case 

In the previous section, we discussed models with a 
fixed physical dissipative scale. The purpose of this sec- 
tion is to demonstrate that without a fixed dissipation 
scale, numerical convergence cannot be reached. In other 
words, for large Reynolds numbers, the system can evolve 
qualitatively differently. 

Figurc lTUI shows a sequence in resolution of models with 
zero physical resistivity and viscosity, i.e. models for 
which the numerical dissipation at resolution scale will 
set the Reynolds number. Thus, higher resolution will 
lead to larger Reynolds numbers. For resolution reasons 
we chose the wave number of the interface perturbation 
to be k = 1 in M3 3.2 3.2~T1 Since the condition for fast 
growth of the instability is given by krj w 1, this required 
a larger initial amplitude perturbation and an elongated 
box. Here, we are interested in the (later) turbulent evo- 




Fig. 8. — Logarithmic density maps for models with the field ori- 
ented transversally, i.e. perpendicular to the inflow, at t = 3.8 Myr. 
Left: Three resolutions increasing from top to bottom (N x = 256 
to 1024 by factors of 2) at ca/c s = 0.5. Right: For ca/c s = 1.0. 
The geometric pattern at the boundaries is an artifact caused by 
magnetosonic modes reaching the inflow boundaries. 



lution of the slab, thus we start with k — 4, which allows 
us to reduce the initial amplitude of the perturbation by 
the same factor 4 and therefore considerably extends the 
spatial range in which the slab can develop. 
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Fig. 9. — Logarithmic magnetic energy maps of models with 
the field oriented transversally, i.e. perpendicular to the inflow, 
at t = 3.8 Myr. Left: Three resolutions increasing from top to 
bottom (N x = 256 to 1024 by factors of 2) at c A /c s = 0.5. Right: 
For ca/c 3 = 1.0. 



At our lowest resolution (N x — 256), we essentially get 
a laminar behavior: the pairwise field reversal regions 
are stretched along the width of the slab but persist to 
the end of the simulation (top). At N x = 1024, the main 
magnetic null regions are accompanied by secondary re- 
gions as a result of additional shear flows. The slab 
is thinner. Increasing the resolution further introduces 
more and more substructure in the slab, especially at the 
"heads" of the slab's perturbations: here, field reversals 
seem to accumulate, leading to additional field dissipa- 
tion. Thus, with higher Reynolds numbers, the slab gets 
more turbulent, and reconnection proceeds not only in 
the two main magnetic null regions, but all throughout 
the slab in small regions. Consequently, turbulent field 
structures inside the slab are dissipated faster, leading to 
a deficit in magnetic pressure inside the slab, and thus 
to a thinner slab. 

While this effect is certainly interesting to note, the sit- 
uation might be less extreme in a truly three-dimensional 
system: a third field component without magnetic null 
could give rise to suffi cient magnetic pressure t o prevent 
reconnection (see also iHeitsch fc Zwe ibcl 200^). In this 
case, small-scale fields entangled by the turbulence in the 
slab could actually lead to additional pressure. 

The pressure profiles fFig. lTl"|) actually tell us a slightly 
more complicated story than that of simple magnetic en- 
ergy dissipation. Pressures were averaged transversally 
(i.e. along the y-axis) and plotted against x, the in- 
flow direction. The magnetic pressure profile stays pretty 
much at a constant level, independent of resolution and 
time. What changes is the kinetic pressure, which drops 
below its inflow value (all panels but bottom). This ef- 
fect gets stronger with increasing resolution. Thus, the 
magnetic field only acts as a dissipation channel for the 
kinetic energy. One could even interpret the kinetic pres- 
sure drop and simultaneous magnetic and internal (red 
lines) pressure rise as an attempt of the system to achieve 
equipartition fFig. I12|) . 

This is not to say that the various pressure compo- 
nents are in equilibrium, as the left panel of Figure 1131 




Fig. 10. — Logarithmic magnetic energy maps for models with 
viscosity v = and resistivity An = 0, i.e. the dissipation scale 
is given by the grid resolution. Resolution increases from top to 
bottom (256 2 to 2048 2 by factors of 2). Since the period in y 
is repeated four times, we show only one quarter of the domain. 
t = 3.8 Myr. 



easily demonstrates. This panel shows the correlation 
coefficient C for the three pairs of pressures, P ma g, Pint, 
and Pkin, for the case without (left) and with (right) 
a physical dissipation scale. Balance between the pres- 
sure components would show up as an anti-correlation, 
whereas a correlation can be interpreted as pressures be- 
ing in phase (and thus driving waves) . Decorrelated pres- 
sures indicate a mixture. Kinetic and magnetic pressure 
decorrelate at higher resolution, pointing to strong recon- 
nection events. Internal and magnetic pressure are only 
slightly correlated (see below), while kinetic and inter- 
nal pressure anti-correlate at late times, because high- 
density regions show more inertia. The models with a 
fixed physical dissipation scale do not show strong reso- 
lution effects. Note, however, that the "dissipation-less" 
models are farther in the dynamical evolution than the 
"controlled" models, because of the different initial con- 
ditions (k y = 4 against k y = 1). 

These results demonstrate that only a fixed dissipa- 
tive scale can guarantee full convergence of the models 
with resolution. Relying on numerical diffusion leads to 
flows with increasing Reynolds numbers as the resolu- 
tion increases and results that depend qualitatively and 
quantitatively on resolution. 

3.3. Growth Rates 

Figure ITTI summarizes the growth of the slab's ampli- 
tude with time. The hydrodynamical growth rates are 
consistent with the analytical predictions (eq. [Q, solid 
straight black line). Saturation sets in when the focal 
points arc shut off from the inflow (see also left col- 
umn of Fig. |SJ). The two highest resolution runs have 
converged also in terms of growth rates (we established 
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Fig. 11. — Transversally averaged pressure profiles for four times 
(in Myr), plotted in ascending time order from bottom to top, 
against x-axis coordinate. Shown are pressure profiles for magnetic 
models with u = and Aq = (see text). Gas is streaming in from 
the left and from the right. See bottom panel for color code and 
line styles. Shown are the total, the kinetic, the magnetic and the 
internal pressures. 




time [Myr] 

Fig. 12. — Average pressure in slab agai nst time. Color coding 
and line style are the same as in Figure ITTI 




Fig. 13. — Correlation of pressures within slab against time. 
Line colors denote correlations, e.g. red stands for C(P m ag, Pkin)- 
Line styles denote resolution. Left: Models without fixed physical 
dissipation scales. Right: Viscous and resistive scales have been 
set. 



detailed convergence in M3 3.2 3.2T|l . Lower resolutions 
lead to slightly smaller amplitudes initially. The low- 
est resolution model is already resolved high enough to 
reproduce the laminar result, but the specified physical 
viscosity is too small to guarantee convergence with re- 
gards to turbulent substructure in the slab. Thus, when 
increasing to the next higher resolution, the unresolved 
turbulence leads to less converged growth rates. Note 
that the growth rate as given in equation Q does not in- 
clude the effect of turbulence generation within the slab, 
alt hough the possi ble effects of turbulence are discussed 
bv lVishniacl (|1994|) . 

The dotted line in Figure [21 and all subsequent figures 
of the same type denotes the growth of the amplitude of 
an unperturbed shock-bounded slab of width A(t), 



A(t) = 2 c s t (M/2 + (1 + (Af/2) 2 ) 1/2 ) 



(4) 



Thus, all models except for those with field strength 
ca/cs = 2 show a faster growth of the slab ampli- 
tude than just the shock-bounded expansion. Obvi- 
ously, equation only gives a rough estimate of the field 
strength required, and in fact underestimates the effect 
of the field. Our crude approximation treats the slab as 
a solid wall, allowing to split the pressure exerted by the 
inflow on the wall along the normal and tangential direc- 
tions without any losses. This is most likely unrealistic, 
i.e. the tangential component will generally be smaller, 
thus requiring a smaller field to balance it. 

With increasing field strength, the ordering of the 
curves with respect to resolution is inverted: now, the 
lowest-resolution runs show the fastest growth. Never- 
theless we get convergence for the two highest resolution 
simulations. This inversion is a consequence of the slight 
increase in the Reynolds number with respect to the runs 
at lower resolution. As discussed above, higher Reynolds 
numbers lead to more turbulence and more field rever- 
sals. Thus, the field is dissipated faster, leading to a 
pressure deficit in the slab. 

A vertical field leads to a strongly reduced growth rate 
or suppresses the instability completely, depending on 
its strength (Figure EJ. The NTSI arises, however, for 
weaker field strengths, and the bending mode in the slab 
persists, thus allowing material to be funneled towards 
the focal points. In the presence of self-gravity or a strong 
thermal instability (e.g. that provided by atomic line 
cooling), these density enhancements could then frag- 
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Fig. 14. — Slab amplitude against time for four models at zero 
field strength ca/c s = to highest field strength Cj\/c s = 2. Line 
styles stand for resolution. The st raight solid black line denotes the 
analytical solution (eq.0, Vishniac 1994). The dotted line denotes 
the expansion of a shock-bounded unperturbed slab (eq. HT). 
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Fig. 15. — Slab amplitude against time for models with trans- 
verse fields. Thin lines denote the corresponding models with fields 
oriented along the inflow for comparison. Again, t he straight soli d 
black line denotes the analytical solution (cq. 0, Vishniac 1994). 
The dotted line denotes the expansion of a shock-bounded unper- 
turbed slab (eq. HIV 



ment and generate further substructure, despite the ini- 
tially unfavorable field orientation. 

The amplitude growth corresponding to Figure I1UI is 
shown in Figure ^j] With higher resolution, saturation 
sets in earlier in the magnetic runs, and the amplitudes 
even decrease with time, mirroring the loss of pressure in- 
side the slab due to (resistive) dissipation. Note that the 
slab amplitude now grows faster than that of an unper- 
turbed shock-bounded slab (see dotted line in Fig. I16fl 
only at times t < 1 Myr. As Figure ^| already had 
made us suspect, saturation sets in earlier for models 
with higher k. 

3.4. Field-Density Relation 

The field-density relation B(n) is often used as an 
observational measure for the dynamical importance of 
magnetic fields in the interstellar medium. It describes 
the mass-loading of field lines, and is related to the 



12 3 

time [Myr] 

Fig. 16. — Slab amplitude against time for models with initial 
perturbation wave number k = 4 and without fixed resistive scale. 
Higher resolution leads to saturation of the growth rate at lower 
amplitudes. Th e straight solid black line denotes the analytical 
solution (cq. 0, Vishniac 1994). The dotted line denotes the ex- 
pansion of a shock-bounded unperturbed slab (cq. HT), 



mass-to-flux ra tio (e.g. iMouschovias &: Spitzerl Il976t 
iShu et aLlll987|) quantifying the importance of magnetic 
fields in a gravitationally dominated cloud. Microscopi- 
cally, the mass loading of field lines can be changed in two 
ways: by Ohmic diffusion and by ion-neutral (or ambipo- 
lar) drift. However, the two parameters controlling the 
degree to which the interstellar magnetic field is frozen to 
the gas are huge. The magnetic Reynolds number KeM 
is the ratio of the Ohmic diffusion time to the dynamical 
time, and is typically of order 10 15 — 10 21 . The sec- 
ond parameter, the ambipolar Reynolds number Kbadi 
is the ratio of the ion-neutral drift time to the dynamical 
time. This number is typically many orders of magnitude 
less than the first one, and can approach unity in dense 
molecular gas. These numbers suggest that the mag- 
netic field should be nearly perfectly frozen to the plasma 
component of the gas, and generally well frozen to the 
neutrals, except in the densest, nearly absolutely neutral 
regions. Thus, the B(n) relation is determined primarily 
by dynamical rather than by microscopic processes. Pa- 
rameterizing the relation by q = dhiB/dhin, one finds 
q = 1 for compression perpendicular to B, q = 2/3 
for isotropic compression, q = 1/2 for self-gravitating, 
magn etically sub-critical clouds ijFiedler & Mouschoviasl 
1993), and q = for compression parallel to B. 

Observations of the B(n) relation in molecular gas 
indeed show that the stro ngest fields are associated 
with the densest gas (e.g. ICrutcherlll999|) . However, 
in the more diffuse ISM, the B{n) relation is consis- 
tent with q m over 3 orders of ma gnitu de in n 
ijTrolann 1 fc Heilesl ll 986. He iles Rr. Trol^ndlEinl . Thus, 
processes beyond microscopic diffusion are required to 
decouple field and density. The possibility of accelerat- 
ing the dec oupling th rough turb ulent transport has been 
expl ored bvlZweibell ll2002l) andlHeitsch et all J2004I) (see 
alsolKim fc Diamondll2002llFatuzzo fc Adamsll2002l and 
ILi &: Nakamural |2004j) . Numerical evidence for a weak 
B(n) relation includes [P adoan fc Nordlundl l|1999() and 
Ide Avillez fc BreitschwerdtJ l|2005|) . 

Figure El shows scatter plots of log B against log n 
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for four models, each measured at t = 4 Myr. The top 
row shows models with the field parallel to the inflow 
(denoted by ca x in the label), while the field is oriented 
transversally, or perpendicularly to the inflow in the bot- 
tom row (denoted by ca v )- Idealized scalings are indi- 
cated by the dashed lines. 

The most striking difference between the top and the 
bottom row is that for the configuration where the field 
is parallel to the inflow, field and density seem to corre- 
late more strongly than for the configuration where the 
field is perpendicular to the inflow. This might seem sur- 
prising at first. After all, one would expect the field to 
be correlated least with density if the gas is compressed 
along the field lines (top row), and to be correlated most 
with density when the gas is compressed perpendicularly 
to the field lines (bottom row). However, looking back at 
Figures [7| and El on the one hand, and Figure EI] on the 
other, we realize that the models with fields parallel to 
the inflow generally develop turbulence, leading to field 
line stretching. Thus, in the (denser) slab, the field gen- 
erally will be stronger. Why, then, is there close to no 
correlation observable in the bottom row of Figure I17f 
The answer is hidden in the way we set up the initial con- 
ditions. To avoid the generation of strong MHD waves, 
we kept the field uniform, despite the fact that the flow 
collision interface is strongly perturbed (see top row of 
Fig- El for the initial conditions). Thus, when gas is de- 
flected at the flanks (point "0" in Figure EJ), it is essen- 
tially free to move along the field lines, i.e. transversally, 
but is increasingly prevented from continuing its trip to- 
wards the slab, because the magnetic pressure is increas- 
ing (note that the bulk magnetic field strength is higher 
in the bottom row of Fig. ll7J) . Thus, the density can take 
on any values in the slab, while the field value is given 
by the ram pressure. A weaker field (lower right panel vs 
lower left panel of Figure ITTjl reduces this effect, allowing 
some scatter in the field strength, and, indeed, checking 
Figure El for this case, the slab actually shows substruc- 
ture and is allowed to bend. Note that there is a strict 
d In B/d In n = 1 scaling for the strong field model (lower 
left panel), which results from the initial compression of 
the field lines. 

While it is hard to generalize the results of our two- 
dimensional models, they demonstrate that the B(n) re- 
lation is strongly influenced by the geometry of the fields 
and the gas flows. Even compression perpendicular to 
the field lines can lead to a nearly complete decoupling 
of field and density - as long as the gas flow is given the 
chance to break the symmetry. 

In a sense, we expect the "geometrical" mechanism 
decoupling the field from the density (lower row of Fig- 
ure I17|) to compete with the turbulent transport of the 
magnetic field: both act on dynamical (flow) timescales. 
For fields oriented perpendicularly to the inflow, our 
models do not develop any substantial turbulence (this 
might also be due to a too small Reynolds number). On 
the other hand, the models with field parallel to the in- 
flow direction do generate some turbulence. For these 
models, B and n decorrelate only at higher density val- 
ues, corresponding to small scales, while the lower den- 
sity regions show a reasonably well established correla- 
tion between B and n. 

4. SUMMARY 



Th e Non-linear Thin Shell Instability (NTSI. iVishniacl 
Il994|) is expected to occur in expanding shells, shocks or 
colliding gas streams. Previous studies have addressed 
the evolution of the NTSI under hydrodynamical condi- 
tions, including gravity and cooling. We have presented 
a numerical study of the NTSI including magnetic fields. 
We have established that our numerical method is well 
suited to tackle the problem. We have found that the ef- 
fects of magnetic fields on the NTSI can be summarized 
as follows: 

(1) Fields principally tend to weaken or even suppress 
the NTSI. We further distinguish between two cases: (i) 
fields aligned with the inflow resist the transverse mo- 
mentum transport - which is the main driving agent of 
the NTSI - via the magnetic tension force; (ii) fields per- 
pendicular to the inflow lead to a stiffer equation of state. 
If ca ~ u, the NTSI is suppressed. However, even for 
transverse fields, substructures can form within the slab, 
which can serve as fragmentation seeds in the presence 
of thermal instabilities or self-gravity 

(2) A fixed physical scale both for viscous and resistive 
dissipation is necessary to reach numerical convergence. 
When relying on numerical dissipation at the resolution 
scale, the Reynolds number will increase with resolution, 
leading to a more turbulent environment and thus to 
results which qualitatively and quantitatively depend on 
resolution CFigsH ITO! and IT1)|) . 

(3) At larger Reynolds numbers, turbulent rcconnec- 
tion plays a role in the turbulent dense slab generated 
by the NTSI. Magnetic energy is therefore dissipated at 
higher rates, leading to a pressure deficit in the dense 
slab. The magnetic field acts as a dissipation channel 
(Figs ITTl and IT5jl . 

(4) Although the energies (or average pressures) seem 
to show a tendency of the system to evolve towards 
equipartition, pressures do not balance locally within 
the slab fFigs. Upland 1 13[l . Correlated pressures lead to 
waves, i.e. the slab's inner structure is highly dynamical. 

(5) The relation between field and density is, at best, 
weak in all models fFig. 1171) . Models with fields paral- 
lel to the inflow exhibit a stronger B(n) correlation than 
models with fields oriented perpendicularly to the inflow. 
The main reason for this is the generation of turbulence, 
which leads to field line stretching and thus field amplifi- 
cation within the denser slab. Fields oriented perpendic- 
ularly to the inflow allow instreaming material to move 
laterally, permitting the field and density to decorrelate. 

Our isothermal models only allow a limited exploration 
of the effect of fields on colliding flows in a thermally or 
gravitationally unstable medium. Clearly, substructure 
can form in the slab under most conditions, providing 
potential seeds for thermal or gravitational instabilities. 
Thus, to establish the role of magnetic fields for molec- 
ular cloud formation in the colliding flow scenario, the 
thermal and gravitational effects have to be addressed. 
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